##=============================================================================
## Appendix: Table 7
##=============================================================================

##-----------------
# clear environment
rm(list=ls())
options(stringsAsFactors = FALSE, scipen = 999)
# source("R/functions.R")

##--------
# Set seed
seed <- sample.int(.Machine$integer.max, 1)
set.seed(seed)

##-------------
# Load Packages
ipak <- function(pkg){new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}

packages <- c("tidyverse", "lme4", "sjPlot", "ordinal", "texreg")


ipak(packages)

##---------
# Load Data
#setwd("~/Dropbox/Ukraine2022WartimeSurvey/Paper_peace/final_version_oct_2024/replication-scripts")
load("clean_ukraine_data.RData")

##-----
# Fit 1
dat <- dat %>% 
  mutate(treat1 = ifelse(treated == 1, 1, 0)) %>% 
  mutate(treat2 = ifelse(treated == 2, 1, 0))

dat$save_lives <- factor(dat$save_lives)

fit1 <- clmm(save_lives ~ 1 + wave + treat1 + treat2 + survival_mindset + per_idx_invcov + 
               fam_idx_invcov + frontline_b + displaced + sex + rus_ethnic +
               age + edu + job + ln_viina_two_months + finances +
               (1 + wave + survival_mindset + per_idx_invcov + fam_idx_invcov | 
                  adm1_en), data = dat)

sjPlot::tab_model(fit1)
#performance::icc(fit1)
#sjPlot::plot_model(fit1)
#AIC(logLik(fit1))

##-----
# Fit 2
fit2 <- glmer(compromise ~ 1 + wave + treat1 + treat2 + survival_mindset + per_idx_invcov + 
                fam_idx_invcov + frontline_b + ln_viina_two_months +
                displaced + sex + rus_ethnic + age + edu + job + finances + 
                (1 + wave + survival_mindset + per_idx_invcov + fam_idx_invcov | 
                   adm1_en), 
              data = dat, 
              family = binomial,
              control = glmerControl(optimizer = "bobyqa"))

sjPlot::tab_model(fit2)
#sjPlot::plot_model(fit2)
#performance::icc(fit2)
#AIC(logLik(fit2))

##-----
# Fit 3
fit3 <- lmer(compromise_idx_terr_z ~ 1 + wave + treat1 + treat2 + survival_mindset + 
                per_idx_invcov + fam_idx_invcov + frontline_b + 
                ln_viina_two_months + displaced + rus_ethnic + sex + age + 
                edu + job + ln_viina_two_months + finances + 
                (1 + wave + survival_mindset + per_idx_invcov + 
                              fam_idx_invcov | adm1_en), data = dat)

sjPlot::tab_model(fit3)
#sjPlot::plot_model(fit3)

##-----
# Fit 4
fit4 <- glmer(ukr_no_nato ~ 1 + wave + treat1 + treat2 + survival_mindset + per_idx_invcov + 
                fam_idx_invcov + frontline_b + ln_viina_two_months +
                displaced + sex + rus_ethnic + age + edu + job + finances + 
                (1 + wave + survival_mindset + per_idx_invcov + fam_idx_invcov | 
                   adm1_en), 
              data = dat, 
              family = binomial,
              control = glmerControl(optimizer = "bobyqa"))

sjPlot::tab_model(fit4)
#sjPlot::plot_model(fit4)

##-----
# Fit 5
fit5 <- glmer(rus_lang_schools ~ 1 + wave + treat1 + treat2 + survival_mindset + per_idx_invcov + 
                fam_idx_invcov + frontline_b + ln_viina_two_months +
                displaced + sex + rus_ethnic + age + edu + job + finances + 
                (1 + wave + survival_mindset + per_idx_invcov + fam_idx_invcov | 
                   adm1_en), 
              data = dat, 
              family = binomial,
              control = glmerControl(optimizer = "bobyqa"))

sjPlot::tab_model(fit5)
#sjPlot::plot_model(fit5)

#screenreg(list(fit1, fit2, fit3, fit4, fit5))

##----------------------
# Export models to LaTeX
coef_names <- c("Survey Wave",
          "Treatment 1",
          "Treatment 2",
          "Survival Mindset",
          "Personal Victimization",
          "Family Victimization",
          "Frontline Oblast",
          "Displaced",
          "Sex",
          "Ethnic Russian",
          "Age",
          "Education",
          "Employment",
          "Ln(VIINA)",
          "Socio-Economic Status",
          "1|2",
          "2|3",
          "3|4",
          "(Intercept)")

mod_names <- c("Save Lives (ordinal)", 
                       "Compromise (binary)",
                       "Compromise (terr. index)",
                       "Compromise (UKR no NATO)",
                       "Compromise (RUS lang schools)")

note <- "Insert footnote here."

texreg(list(fit1, fit2, fit3, fit4, fit5), 
       dcolumn = TRUE,
       omit.coef = "\\(Intercept\\)",
       custom.coef.names = coef_names,
       custom.model.names = mod_names,
       custom.note = str_c("Note: ", note),
       caption.above = TRUE,
       caption = "Regressions of Occupational Prestige",
       digits = 3, 
       leading.zero = TRUE,
       fontsize = "small")

rm(list = ls())
##=============================================================================
## End of File
##=============================================================================